Analytic Design Technique for 2D FIR Circular Filter Banks and Their Efficient Implementation Using Polyphase Approach

This paper proposes an analytical design procedure for 2D FIR circular filter banks and also a novel, computationally efficient implementation of the designed filter bank based on a polyphase structure and a block filtering approach. The component filters of the bank are designed in the frequency domain using a specific frequency transformation applied to a low-pass, band-pass and high-pass 1D prototype with a specified Gaussian shape and imposed specifications (peak frequency, bandwidth). The 1D prototype filter frequency response is derived in a closed form as a trigonometric polynomial with a specified order using Fourier series, and then it is factored. Since the design starts from a 1D prototype with a factored transfer function, the frequency response of the designed 2D filter bank components also results directly in a factored form. The designed filters have an accurate shape, with negligible distortions at a relatively low order. We present the design of two types of circular filter banks: uniform and non-uniform (dyadic). An example of image analysis with the uniform filter bank is also provided, showing that the original image can be accurately reconstructed from the sub-band images. The proposed implementation is presented for a simpler case, namely for a smaller size of the filter kernel and of the input image. Using the polyphase and block filtering approach, a convenient implementation at the system level is obtained for the designed 2D FIR filter, with a relatively low computational complexity.


Introduction
The technology and architecture of modern image sensors and sensing techniques have evolved dramatically in recent years, driven by the ever-demanding requirements and challenges of this field.For instance, aerial or satellite image sensors for remote sensing must provide clear and low-noise images, with high spatial resolution, either in visible, infrared or microwave domains.In order to provide accurate and relevant information, images acquired by sensors have to be pre-processed using various restoration and enhancement techniques.Various digital filters and filter banks may be used in image analysis and feature extraction tasks, for instance, to decompose the image into several subband components in order to extract relevant details, etc.These are also useful in the automotive field, for rapid feature extraction in real-time computer vision applications, for instance, in driver assistance systems and autonomous driving vehicles.
Along with the unprecedented development of the digital signal processing field, 2D filters have been thoroughly investigated by many researchers, owing to their essential applications in image processing, and various techniques for their design have been elaborated [1].Analytical design methods rely on 1D prototypes with specified shapes and parameters; applying various frequency transformations, they lead directly to the desired Sensors 2023, 23, 9851 2 of 21 2D filters.The major advantage of the analytical approach is that a closed-form frequency response is derived, and the 2D filter results are parametric and therefore adjustable.
A large variety of 2D filters, both of FIR and IIR type, with various characteristics and shapes have been developed, with each type of filter having specific applications in the image processing field.One of the best-known methods, widely used in the design of 2D FIR filters with various shapes is the McClellan transform [2,3].More recent papers approaching computationally efficient 2D FIR filter design techniques based on frequency transformations are [4,5].The efficient, low-complexity design of 2D FIR filters and the implementation using Farrow structure are described in papers like [6,7].Other relevant recent papers on efficient 2D filter design are [8][9][10].
Filters with circular-shaped frequency response have also been widely used owing to their capabilities in image analysis; various design techniques have been proposed for circular filters (CF) in early papers such as [11][12][13].Circular filters find applications in texture segmentation and classification [14].A recent advanced application of circular Gabor filters in SAR interferograms is described in [15].Two-dimensional filter banks of various types were extensively used in important applications, like texture segmentation and classification or various feature extraction tasks.Such filter banks decompose the frequency spectrum of the image into a number of sub-bands.Two-dimensional filter banks are widely used in fundamental applications such as sub-band coding and compression of images and video sequences.Separable 2D filter banks are obtained by cascading 1D filter banks, and data are processed in each dimension separately.Compared to separable filters, filter banks with nonseparable 2D filters are more flexible and versatile, offering superior performance for imposed specifications.However, their design is substantially more difficult than for separable filter banks [16].As detailed in the comprehensive review [16], the 2D filter banks currently used are mainly directional, with specific shapes in the frequency plane, such as square (diamond), parallelogram, wedge/fan filters, etc. Multidimensional stable, perfect reconstruction filter banks are also developed in [17].
Directional filter banks (DFBs) with an arbitrary number of sub-bands [18] or arbitrary frequency partitioning [19] have been proposed.A class of multiresolution DFBs is developed in [20].Multidimensional DFB, multiscale pyramids and the surfacelet transform were introduced in [21].A very recent application of DFBs was proposed in [22], namely fingerprint image quality assessment.The fingerprint image is decomposed into subbands using the DFB, and similarity between the different subbands is used to calculate the fingerprint image quality.Regarding methods to reduce computational complexity and increase processing speed, the fast block implementation of 2D digital FIR filters was proposed in early papers such as [23].A high-performance 2D parallel block-filtering system for real-time applications was presented in [24].The steerable pyramid, a well-known multiscale structure for image decomposition was proposed in the early paper [25].More recent papers describe specific applications of other two important multiscale architectures, namely the Laplacian pyramid [26] and the wavelet pyramid [27].
Some very recent works propose advanced algorithms implemented on various convolutional neural networks to solve complex image-processing tasks.For instance, in [28], a novel deep-feature model has been proposed for coastal wetland classification using multisource satellite remote sensing data.In [29], multi-scale features from coarse-to-fine receptive field level are extracted, with applications in super-resolution.An advanced algorithm for effective pathology classification from hyperspectral medical images is proposed in [30].A novel multi-focus image fusion method based on sparse representation and local energy is introduced in [31], which uses the shearlet transform to decompose the source images into low-and high-frequency sub-bands.
The first author of this paper has also proposed various analytical design techniques for 2D filters in previous works [32][33][34][35].Directional IIR filters based on Gaussian and wide-band prototypes were designed in [32].A useful application of the directional filters in [32] is the detection of straight lines with specified orientation from images; this feature extraction capability may be useful in the computer vision field.Adjustable, parametric 2D digital IIR filters with elliptical and circular symmetry are proposed in [33].Two versions of circular IIR filter banks and their applications have been described in [34,35].An efficient 2D FIR filter implementation based on a polyphase approach and block filtering is proposed in [36].
In this paper, an analytic design procedure is proposed for a particular class of 2D filter banks, namely 2D FIR Gaussian circular filter banks (CFBs).Two versions of CFBs will be designed, namely a uniform CFB and then a non-uniform (dyadic) CFB, each with a specified number of component filters.As a prototype, a 1D low-pass filter with a Gaussian frequency response and specified selectivity is chosen; its frequency response is easily approximated by a trigonometric polynomial, with an imposed precision, using a simple Fourier series expansion.By a simple shifting to a given peak frequency, the band-pass filters of the FB prototype are also derived.Once the prototype FB is obtained, a 1D to 2D frequency mapping derived from the McClellan transform is applied [2,3], which leads directly to the desired circular filters of the CFB.The non-uniform (dyadic) CFB is designed in a similar manner.The filters' characteristics result in an accurate circular shape, with some distortions near the frequency plane margins.Next, as an application example, a grayscale test image is applied to the CFB, obtaining a set of subband images.Summing back all these images, the original input image is reconstructed almost perfectly, which suggests a potential use in an alternative subband coding scheme.
A novel, efficient implementation solution is also proposed for the 2D FIR filters of the designed CFB, which continues the method from previous work [36].Our implementation uses a polyphase decomposition of a given 2D filtering operation with large kernel size and a block filtering with smaller size matrices.
The paper is organized as follows: Section 2 presents the proposed analytical design procedure, first deriving the uniform and non-uniform prototype FB, then applying the frequency mapping and obtaining the frequency responses of the 2D CFBs.In Section 3, an example of image analysis is given using CFB by decomposing it into subband images.The novel implementation technique based on the polyphase and block filtering approach is described in Section 4. Discussions regarding the computational complexity of the proposed implementation are included in Section 5. Finally, conclusions are drawn in the last section.

Analytical Design Technique for 2D Circular FIR Filter Banks
A novel analytical design procedure is proposed for a class of 2D FIR circular filters.This design technique starts from an imposed prototype with specified parameters (peak frequency, bandwidth), to which a 1D to 2D frequency transformation is applied, leading to the desired 2D filters.In order to obtain through frequency transformation, the desired 2D circular filter bank, first a 1D prototype filter bank must be derived.A Gaussian-shaped filter was chosen as prototype, due to its useful property of scalability on the frequency axis.

Approximation of the Gaussian FIR Filter Prototype Using Fourier Series
The Gaussian filter in the frequency domain has the well-known expression G(ω) = exp −σ 2 ω 2 /2 , where σ is the dispersion parameter; for a simpler form, easier to handle, the substitution p = σ 2 /2, or equivalently σ = 2p, will be used.Thus the Gaussian low-pass filter function takes the more convenient form G LP (ω) = exp −p • ω 2 , where p will be referred to as selectivity or scaling parameter.Considering a periodic function with period 2π and regarding the LP Gaussian function as a generating pulse, the following expression H LP (ω) will be easily obtained, which is the Fourier series expansion of the Gaussian G LP (ω) up to a given order N: • cos nω = H LP (ω) (1)     Sensors 2023, 23, 9851

of 21
From this Gaussian LP prototype, a band-pass (BP) prototype is easily produced by shifting the Gaussian laterally around the frequencies ±ω 0 : Directly using Expressions (1) and ( 2) implemented in a Matlab routine, in the following section, the low-pass, band-pass and high-pass components of the desired FIR filter bank prototype are calculated.

Design of a Gaussian Uniform FIR Filter Bank Prototype
Next, a uniform filter bank prototype with 11 Gaussian components will be designed, namely one low-pass filter, nine band-pass filters and one high-pass filter.In this uniform FB, the peak frequencies are equally spaced on the frequency axis.A bandwidth is imposed for the nine band-pass components equal to B = π/10 = 0.1π, while the low-pass and high-pass filters will have each half of this bandwidth, namely B/2 = π/20 = 0.05π.The k-th ideal Gaussian BP filter is produced by shifting the LP prototype to the frequency ω 0,k = k • ω 0 , and will have the following expression: At this point, the scaling parameter p for the imposed bandwidth needs to be calculated.In our case, the filter bandwidth is considered defined at 0.5 of the peak value (at 6 dB).Thus, the characteristics of any two adjacent filters will marginally overlap and will intersect at the value 0.5.Referring to the LP filter G LP (ω) = exp −p • ω 2 , the condition G LP (B/2) = exp −p • B 2 /4 = 0.5 is imposed, otherwise written exp p • B 2 /4 = 2, from which the value for the scaling parameter p is obtained as p = 4 ln 2/B 2 ; since for our filter bank a bandwidth B = π/10 was imposed, the value p = 400 ln 2/π 2 ∼ = 28.1 will be produced.The ideal uniform Gaussian filter bank is plotted in Figure 1a.
  247.118 (cos +0.99492)(cos +0.95454)(cos +0.87546)(cos +0.76089)(cos +0.61554) Finally, the highest component of the FB is the high-pass (HP) filter, which formally has the peak frequency 0    :   247.118 (cos +0.93698)(cos +0.90811)(cos +0.80475)(cos +0.67403)(cos +0.51542) It can be observed that the component filters of the prototype FB whose central frequencies are symmetric with respect to the middle value 2    have symmetric zeros, as is well-known from filter theory.Therefore, the zeros of the HP filter are the zeros of the LP filter with a changed sign; the zeros of the 9th BP filter are the zeros of the first BP filter with a changed sign, etc.Since there is an odd number of filters, the middle filter, namely the 5-th BP filter, with central frequency 0 2 Generally, the k-th band-pass component of the 1D filter bank can be expressed as the following product of first-order factors (where N is the filter order): The uniform Gaussian filter bank designed above is plotted in Figure 1b and it looks very similar to its ideal counterpart in (a), with a low level of ripple.The filter selectivity is given by the scaling parameter value calculated before, namely p = 28.1, with the Fourier series truncated at a number of terms N = 15.The larger the number of terms taken into account, the smaller will be the distortions (ripple, etc.), but the filter matrices will be larger in size and will increase the implementation complexity.
Following the above design procedure, once specifying the desired number of filters of the FB and their peak frequencies, using Equations ( 1) and ( 2), the frequency responses of all the FB components are calculated.As an example, in our case of a uniform FB with 11 components, the frequency responses of a few filters of the 1D prototype filter bank are given below, in factored expression.First, the frequency response of an LP prototype expressed as a truncated Fourier series using (1) has the form: which using trigonometric identities can be further expressed as: As an example, the frequency responses of the first and last BP filter components of the bank are given below, the intermediate BP filters having similar forms: Finally, the highest component of the FB is the high-pass (HP) filter, which formally has the peak frequency ω 0 = π: H HP (ω) = −247.118• (cos ω+0 .93698)(cosω+0 .90811)(cosω+0 .80475)(cosω+0 .67403)(cosω+0 .51542) •(cos ω+0 .33537)(cosω+0 .14123)(cosω − 0.05906)(cos ω − 0.25729)(cos ω − 0. 44535) •(cos ω − 0. 61554)(cos ω − 0. 76089)(cos ω − 0. 87546)(cos ω − 0.95455)(cos ω − 0.99491) (8) It can be observed that the component filters of the prototype FB whose central frequencies are symmetric with respect to the middle value ω = π/2 have symmetric zeros, as is well-known from filter theory.Therefore, the zeros of the HP filter are the zeros of the LP filter with a changed sign; the zeros of the 9th BP filter are the zeros of the first BP filter with a changed sign, etc.Since there is an odd number of filters, the middle filter, namely the 5-th BP filter, with central frequency ω 0 = π/2, has no pair, and its transfer function, as expected, has pairs of complementary zeros: Generally, the k-th band-pass component of the 1D filter bank can be expressed as the following product of first-order factors (where N is the filter order): The uniform Gaussian filter bank designed above is plotted in Figure 1b and it looks very similar to its ideal counterpart in (a), with a low level of ripple.

Design of a Gaussian Non-Uniform FIR Filter Bank Prototype
In image analysis, mainly in multirate signal processing, non-uniform filter banks are also currently used.Next, using the method described in Section 2.1 a non-uniform, more specifically a so-called dyadic filter bank will be designed.Such an FB has the property that the bandwidths of the component filters increase proportionally to their peak frequencies, such that generally the ratio between bandwidth and peak frequency remains constant; these filters are also known as constant-Q filter banks.
As a further remark, in previous papers [32][33][34][35], various 2D filters were designed using another efficient procedure, namely the Chebyshev series, which has the advantage of yielding a uniform and efficient approximation for a given function, with equal error along the whole specified range of values.The symbolic calculations are performed in the MAPLE software (version MAPLE 2018), and a change of frequency variable is first required, before effectively deriving the approximation.However, the major drawback of this method is that it is not parametric; it does not have a closed form as in the case of the Fourier series method, therefore is more laborious; for each specified value of selectivity parameter p, the calculation must be carried out in a symbolic calculation software.Therefore, in this paper, the Fourier series approximation was preferred.
    The characteristics of this non-uniform FB are plotted in Figure 2.

Gaussian Circular FIR Filter Bank Obtained Using Frequency Transformation
Once specified a convenient 1D prototype with the frequency response H p (ω), a 2D circular filter H(ω 1 , ω 2 ) is produced by applying to the given prototype the 1D to 2D frequency transformation ω → ω 2 1 + ω 2 2 : The function cos ω 2 1 + ω and can be approximated by the following expression, which is a simple particular case of the McClellan transform, currently used in 2D FIR filter design [2,3,36]: The Expression (18) is in fact the discrete space Fourier transform (DSFT) of the matrix C. Next, a zero-phase FIR filter H P (ω) is considered, whose frequency response is given by the trigonometric polynomial expression [36]: At this point, the trigonometric identities for cos kω (k = 1 . . .R) can be used, and thus the following polynomial expression is produced in powers of cos ω [36]: where ( 19), (20) b 0 , b k , c 0 , c k are polynomial coefficients.Applying frequency mapping (18), the frequency response of the 2D circular filter will become [36]: where C(ω 1 , ω 2 ) = cos ω 2 1 + ω 2 2 as given in (18).Therefore, by a straightforward substitution of cos ω by the circular cosine function in the prototype H P (ω), the 2D filter frequency response is produced directly.Next, supposing that the frequency response H P (ω) is decomposed into first-order and second-order factors in variable cos ω, and achieving the above substitution in all factors of H P (ω), the circular filter frequency response H(ω 1 , ω 2 ) is finally derived in factored form: where C is a concise notation for the two-variable function C(ω 1 , ω 2 ) and k is the constant resulting from factorization.Since the specified prototype is expressed as a product of elementary factors, the circular filter frequency response will also become directly factored, which is an essential advantage in actual implementation.Thus, the large kernel H corresponding to H(ω 1 , ω 2 ) can be expressed simply as a discrete convolution of small matrices (of size 3 × 3 or 5 × 5): The matrix expression ( 23) is related to the factored frequency response (22).Using the 3 × 3 matrix C in (17) and considering also (22), each of the matrices C i of size 3 × 3 in (23) is derived by adding coefficient b i , which appears in the first-order factors in (22), to the center element in matrix C. Thus, the matrix D j (5 × 5) becomes: where C 0 is a null matrix of size 5 × 5 with central element of value one; C 1 (5 × 5) is produced by the boarding matrix C (size 3 × 3) with zeros; here the symbol * denotes convolution.Thus, the frequency response of each CFB component is directly derived by substitution.Correspondingly, the overall kernel matrix H of the filter will be given by an expression similar to (23), but in our particular case with only first-order factors, as in (10), it becomes: The filters of the designed 1D prototype filter bank are of order 15; it follows that the corresponding 2D circular filters derived through the above transformation have kernel matrices relatively large, of size 31 × 31.Such a large matrix will be implemented efficiently using a polyphase approach described in Section 4.
As a remark, all the component filters of the designed FB are non-separable, except the LP filter.Indeed, it is easy to see that the circular LP Gaussian filter is separable as a product of two Gaussian LP filters on the two frequency axes: A very important advantage of the proposed FB is that the filters' transfer functions are real-valued (zero-phase), therefore they will not introduce any phase distortions; this will be visible in the simulation results given in the following section.
The 1D prototypes, frequency characteristics and corresponding contour plots for all 11 filters of the circular filter bank are displayed in Figures 3 and 4. It is easily observed that up to the 6th band-pass filter, the characteristics are visually almost perfectly circular.For the higher band-pass filters, the characteristics have a more pronounced deviation from circularity, tending to the shape of a rounded square.The filter with the highest frequency (ω 0 = π) has almost a square shape.This effect of distortion from circularity is well known when applying the frequency mapping (18), the simplest form of the McClellan transform, and could be corrected only by using a more accurate approximation of the circular cosine; however, this would imply a higher complexity of the filters (larger kernel matrices) and a more difficult implementation.The characteristics and contour plots of the five component filters of the non-uniform (dyadic) CFB derived from the 1D prototype filters designed in Section 2.3, with frequency responses given by ( 11)-( 15) are displayed in Figure 5, and it can be observed that they have a good circular symmetry.The characteristics and contour plots of the five component filters of the non-uniform (dyadic) CFB derived from the 1D prototype filters designed in Section 2.3, with frequency responses given by ( 11)-( 15) are displayed in Figure 5, and it can be observed that they have a good circular symmetry.The characteristics and contour plots of the five component filters of the non-uniform (dyadic) CFB derived from the 1D prototype filters designed in Section 2.3, with frequency responses given by ( 11)-( 15) are displayed in Figure 5, and it can be observed that they have a good circular symmetry.

Image Analysis Using the Designed Circular Filter Banks
In this section, examples of image analysis using the uniform and dyadic CFBs designed before are presented.First, the grayscale test image in Figure 6a is considered, of size 399 × 399 pixels, representing a group of trees without foliage; this image was chosen as it has a lot of fine details, represented by the tree ramifications into thinner and thinner twigs.This image is filtered by applying all the 11 components of the designed uniform CFB (one LP filter, nine BP filters, one HP filter); it can be considered that our test image is decomposed into sub-bands using the analysis CFB designed before.
images, as shown in Figure 8f-h.Summing up all the sub-band images leads to an image very similar to the original one.The original image is displayed in Figure 6a.The image obtained at the output of the narrow LP filter is (b), and it can be observed that it is very blurred, the fine details (thin twigs) are no longer visible.The images obtained from the first five BP filters are shown in (c-g), respectively, and contain details corresponding to the selected bandwidth.The image (h) is produced at the output of the HP filter and contains the highest frequencies, corresponding to the finest details.
The original image was converted into "double" format and its pixel values were rescaled to the range [0, 1] for MATLAB processing.The image produced at the output of LPF has the overall mean pixel value 0.529; for all the other 10 images (produced at the outputs of BP filters and HP filter), the mean pixel value is very close to zero, as expected, since these filters eliminate the zero-frequency component corresponding to mean value.In Figure 6i-l, it is shown how the original image is reconstructed by adding the component images into which it was decomposed.Thus, image (i) is produced by adding the first two components (LP and BP1); image (j) is produced as a sum of the first three components (LP, BP1, BP2); image (k) is produced by adding the component BP3.Finally, by summing all the 11 components, the image (l) is produced, which visually is very similar to the original image, showing all the fine details very clearly.
These simulations prove that the designed CFBs (uniform and non-uniform) could be practically used as analysis filter banks for decomposing a given image into sub-band images.However, the rigorous mathematical conditions required will have to be further investigated in future work.
The energy of each component sub-band image can also be evaluated using the well- , where the image is of size M × N and p ij is the current pixel value; the expression of the relative energy can also be given as a percentage: Calculating the energies of the 11 filtered images resulting at the output of the designed CFB, the values given in Table 1 are easily found; summing these values, it can be verified that they add up to approximately 1 in normal values, or 100% in percentages.It can be observed that almost 56% of the image energy is contained in the low-pass component (in the frequency domain around zero, with radius 0.1π), while almost 85% is contained in the first four components (within a 0.4π radius), at the output of LP filter and first three BP filters.The relative energies of the sub-band images decrease almost uniformly; as an exception, ER9 > ER8, and ER11 > ER8, ER9, ER10.The highest frequencies in the image give less than 2% of the total image energy.These relative energy values are summarized in Table 1 and represented graphically in the chart from Figure 7.As a remark, the designed circular filter banks are rotation invariant; the image spectrum is separated into concentric, ring-shaped regions, with frequencies increasing while image energy is generally decreasing, from the center to the margins of the frequency plane.Due to rotational invariance, the decomposition coefficient and energy in each subband remain more or less constant.This property is very useful in specific feature extraction and classification tasks in image processing.
A similar experiment was performed using the dyadic circular filter bank with five components shown in Figure 5, applied on the same grayscale test image, for comparison.The filtered images obtained at the output of the LP filter, three BP filters and HP filter are displayed in Figure 8a-e.As in the previous example, the original image is then reconstructed by adding the first two and three sub-band images, then all the five sub-band images, as shown in Figure 8f-h.Summing up all the sub-band images leads to an image very similar to the original one.An additional experiment was performed using the same dyadic CFB, applied on another grayscale test image ("Fields"), of size 699 × 699, showing an aerial view of a rural landscape with fields and a river (Figure 9a).The filtered images obtained at the output of the LP filter, three BP filters and HP filter, respectively, are displayed in Figure 9b-f.As in the previous examples, the original image is then reconstructed by adding the first two sub-band images (Figure 9g), then all the five sub-band images, as shown in Figure 9h.Summing up all the five sub-band images yields an image very similar to the original one.Table 2 displays the relative sub-band energies, calculated for both test images, namely "Trees" and "Fields".Again, most of the image energy is contained in the lowest sub-band (corresponding to the LPF); however, the energy distribution clearly depends on the particular image, as was expected, an can be considered a numerical indicator characterizing the sub-band decomposition of a given image.
Regarding the noise suppression issue, the authors did not intend to investigate it in this paper.Of course, noise removal is a very important task in image enhancement and restoration.Considering the nature of the noise (Gaussian, salt-and-pepper, speckle noise, etc.), a specific type of filter should be chosen to remove it optimally.Anyway, An additional experiment was performed using the same dyadic CFB, applied on another grayscale test image ("Fields"), of size 699 × 699, showing an aerial view of a rural landscape with fields and a river (Figure 9a).The filtered images obtained at the output of the LP filter, three BP filters and HP filter, respectively, are displayed in Figure 9b-f.As in the previous examples, the original image is then reconstructed by adding the first two sub-band images (Figure 9g), then all the five sub-band images, as shown in Figure 9h.Summing up all the five sub-band images yields an image very similar to the original one.Table 2 displays the relative sub-band energies, calculated for both test images, namely "Trees" and "Fields".Again, most of the image energy is contained in the lowest sub-band (corresponding to the LPF); however, the energy distribution clearly depends on the particular image, as was expected, an can be considered a numerical indicator characterizing the sub-band decomposition of a given image.
noise removal should be achieved before any further image analysis.For the proposed CFB, since the image spectrum is partitioned into ring-shaped regions corresponding to sub-band images, if the original image was affected by some type of noise, it would be distributed more or less evenly in the sub-band images, mainly in higher frequency bands.Therefore, it would have to be eliminated separately from each sub-band component image, which may be a more difficult task.This issue remains to be studied in future work on this topic.

Polyphase Implementation of the Designed 2D Circular FIR Filters
In the following, a low-complexity implementation is proposed for the 2D FIR circular filter bank previously designed, relying on a polyphase structure of a 2D filtering task with a convolution kernel of relatively large size ( 31 31  ).In order to achieve convolution with such a large kernel, a block processing technique [24,25] and a polyphase decomposition approach will be employed.
As a first step, using sub-expression sharing techniques, a 2D filtering algorithm with a 4 × 4 kernel was elaborated, which is detailed as follows.The kernel of the filter resulting from the design and the input image are decimated by factors 3 and 5, respectively; the polyphase filtering approach is subsequently applied.Using this technique, three output component images are derived, namely 0 Y , given by Equations ( 28)-(30):  Regarding the noise suppression issue, the authors did not intend to investigate it in this paper.Of course, noise removal is a very important task in image enhancement and restoration.Considering the nature of the noise (Gaussian, salt-and-pepper, speckle noise, etc.), a specific type of filter should be chosen to remove it optimally.Anyway, noise removal should be achieved before any further image analysis.For the proposed CFB, since the image spectrum is partitioned into ring-shaped regions corresponding to sub-band images, if the original image was affected by some type of noise, it would be distributed more or less evenly in the sub-band images, mainly in higher frequency bands.Therefore, it would have to be eliminated separately from each sub-band component image, which may be a more difficult task.This issue remains to be studied in future work on this topic.

Polyphase Implementation of the Designed 2D Circular FIR Filters
In the following, a low-complexity implementation is proposed for the 2D FIR circular filter bank previously designed, relying on a polyphase structure of a 2D filtering task with a convolution kernel of relatively large size (31 × 31).In order to achieve convolution with such a large kernel, a block processing technique [24,25] and a polyphase decomposition approach will be employed.
As a first step, using sub-expression sharing techniques, a 2D filtering algorithm with a 4 × 4 kernel was elaborated, which is detailed as follows.The kernel of the filter resulting from the design and the input image are decimated by factors 3 and 5, respectively; the polyphase filtering approach is subsequently applied.Using this technique, three output component images are derived, namely Y 0 , Y 1 , Y 2 , given by Equations ( 28)-( 30): in which the block matrices have the form given below: and O 4×10 , O 10×4 , O 10×7 are zero matrices of size 4 × 10, 10 × 4, and 10 × 7, respectively.Adding the partial results Y 0 , Y 1 and Y 2 given by ( 28), ( 29) and ( 30), the following output vector Y containing 16 samples of the filtered image is obtained: The vector H occurring in Equations ( 28)-( 30) is given below: The main reason for proposing this 2D FIR filtering algorithm was the reduced number of arithmetic operations involved.It is well known that in a direct 2D convolution there is a high degree of redundancy in operations.In a direct 2D convolution there are overlapping blocks of input data; by eliminating these redundant calculations, a significant reduction in the arithmetic complexity will be obtained.The filtering algorithm presented above was produced using a block filtering technique.
At this point, the 2D filtering algorithm discussed above will be extended from an elementary kernel of size 4 × 4 to the case of a 31 × 31 kernel.In order to achieve this and to obtain a parallel implementation, a block processing technique will be used, relying on a polyphase structure.To derive this 2D polyphase structure, a decimation of the kernel matrix with factor 4 will be performed.Before decimation, the kernel was enlarged to have a dimension multiple of 4, in our case 32 × 32, by bordering it with a row and a column of zeros.Decimation by a factor of 5 was also applied to the input image and thus a 25 × 25 input image was produced.
Using a block polyphase decomposition and the previous fast algorithm, the following efficient algorithm was obtained for the computation of the designed 2D FIR filter.given below (where i = 0, 1, 2, 3 and j = 0, 1, 2, 3): For example, the vectors H 00 , H 12 , H 33 generated by the Formula (35) will be: In order to explain the proposed method in an easier way, our demonstration was restricted to a less complex particular situation where the kernel matrix is of size 12 × 12 and the input image is 21 × 21, but the results can be readily extended for the kernel of the circular FIR filter designed above of size 31 × 31, previously extended to size 32 × 32 (by padding with zeros), to be able to achieve the decimation by a factor of 4. The simpler algorithm described above for a 2D filter with a 3 × 3 kernel and 5 × 5 input matrix, can be extended by performing a decimation by factor 4 for the kernel matrix and a decimation by factor 5 for the input matrix.Thus, performing a decimation with factor 4, instead of the kernel of size 12 × 12, 16 matrices of size 3 × 3 are derived.For instance, in the case of H T 01 , applying decimation by 4, the following block matrix of size 3 × 3 will produce: Next, by concatenating the rows of matrix H 01 , the matrix H T 01 is derived from Equation (35).The vector X 2D is also substituted with vector X 2D given by (34).
The vectors X 00 , X 01 , X 02 , X 03 , . .., X 66 , composing the matrix X 2D and related to the input image, are defined through the following general formula: For example, the vectors X 03 , X 31 , X 66 generated by the Formula (38) will be: X 03 = x 14,17 x 14,10 x 14,3 x 7,17 x 7,10 x 7,3 x 0,17 x 0,10 x 0,3 ( i = 0, j = 3) X 31 = x 17,15 x 17,8 x 17,1 x 10,15 x 10,8 x 10,1 x 3,15 x 3,8 x 3,1 ( i = 3, j = 1) X 66 = x 20,20 x 20,13 x 20,6 x 13,20 x 13,13 x 13,6 x 6,20 x 6,13 x 6,6 The vectors X 00 ,. .., X 66 were produced as described below.To explain the idea, our demonstration is restricted to the situation in which the input image is a matrix of size 21 × 21 and a decimation by factor 5 is performed.In doing so, instead of the input matrix of dimension 21 × 21, a number of 77 matrices of size 3 × 3 are obtained.For instance, in the case of X 01 , applying decimation by the factor 7, the following 3 × 3 block matrix is derived: At this point, the rows of matrix X 01 are concatenated, then the resulting vector is reversed and thus the vector X 1,0 is derived from the general Equation (38): X 1,0 = x 15,14 x 15,7 x 15,0 x 8,14 x 8,7 x 8,0 x 1,14 x 1,7 x 1,0 Even if our discussion was restricted to the particular case where the input matrix is 21 × 21 to be easier for the reader to follow our discussion, it is easy to extend it for a more general case.Thus, a 2D FIR filtering operation with a 4 × 4 kernel and a 7 × 7 input matrix was decomposed into 100 1D inner products (FIR filtering operations) using the following equations: Finally, the following output vector is produced: In Equations ( 42)-(44), the block matrices are, respectively: B 0 = A 0 ⊗ U 9 , where the vector U 9 is U 9 = 1 1 1 1 1 1 1 1 1 and B 1 = A 1 ⊗ I 9 , where I 9 is the 9 × 9 identity matrix (with ones on the main diagonal and zeros elsewhere); we also have B 2 = A 2 ⊗ I 9 .The matrices O 4×90 , O 90×36 and O 90×63 are zero matrices of size 4 × 90, 90 × 36 and 90 × 63, respectively.
In order to obtain the above equations, we considered a polyphase decomposition of a 1D filter that can compute four samples in parallel using a decimation factor of 4 as: By extending the Equation (46) to 2D and using sub-expression sharing, we obtained Equations ( 42)-( 44).Although at first sight, the matrix equations describing the proposed polyphase implementation may seem very complex, mainly due to their block structure, they actually lead to a very efficient and economic filtering structure, with a high degree of parallelism and therefore with a low computational complexity in terms of number of arithmetic operations.All these equations were verified in Matlab (version R2017a).

Discussion
The proposed design method for 2D circular filters is entirely analytical, without using any global numerical optimization techniques.Analytical design methods lead to closed-form and parametric filters, with adjustable, tunable frequency responses.To the best of the authors' knowledge, the analytical design of FIR circular filter banks has not been systematically approached previously by other researchers.As a reference to existing works, analytical techniques for designing 2D filters of IIR type with circular frequency response, including CFBs, have been previously proposed by the first author [33][34][35].
The Gaussian filter was chosen as a prototype for the CFB due to its advantages.It is a smooth function that can be easily approximated by a trigonometric polynomial and can be scaled on the frequency axis to adjust its selectivity.For very selective filters, the Gaussian shape is probably the best choice.Its frequency response is zero-phase; since frequency components will not be phase-shifted, image distortions will not occur.The resulting filters have accurate shapes, with negligible distortions.Moreover, they can be approximated efficiently, leading to low-order filters.
The circular filter bank (CFB) designed in our paper can be compared with other types of filter banks, from a qualitative point of view.The comparative discussion will mainly refer to works [25][26][27], as well-known multiscale pyramidal decomposition methods.Our proposed filter banks, like the steerable pyramid [25], have rotation invariance, while the Laplacian pyramid [26] and wavelet pyramid [27] are not rotationally invariant.Another important aspect regards frequency plane partitioning.While the steerable, wavelet and Laplacian pyramids all split the image spectrum into fixed sub-band regions, the proposed circular FB is flexible, in the sense that the bandwidths of the sub-band regions can be chosen wider or narrower, with adjustable selectivity, depending on application.This is due to the scalability of Gaussian-shaped filters along the frequency axis, which allows us to obtain filters with imposed selectivity starting from the same prototype.In Section 3, a uniform CFB with 11 components was generated, partitioning the image spectrum into concentric ring-shaped sub-bands.Moreover, in some applications, the non-uniform (dyadic) filter bank is also useful from the multi-resolution point of view.Since the energy of an image spectrum is mainly contained in the low-frequency region and decreases towards higher frequencies, the dyadic-type CFB allows for a more uniform energy distribution on frequency bands.Also, the proposed polyphase implementation structure for the filter bank has a lower arithmetic complexity than other implementations found in the literature.
A rigorous comparison in terms of performance with other circular filters found in literature is quite difficult to make.Design approaches like circular filters in [11][12][13] are very different from the one proposed here and lead to filters with other characteristics and purposes, so they are quite difficult to compare exactly with our proposed method.
To summarize, the novelty of the proposed CFB consists of an analytic design method (yielding parametric, closed-form expressions of frequency response), frequency scalability, flexible partitioning of spectrum sub-bands, low order due to efficient approximation and low arithmetic complexity due to polyphase implementation.
The proposed novel implementation technique significantly reduces the number of arithmetic operations required.A short comparison can be made between the direct convolution operation and the proposed filtering method in terms of computational complexity.The 2D filtering of an image of size M × N pixels, with an FIR filter with kernel size m × n implies a 2D convolution between a m × n matrix and a M × N matrix.This means that the filter kernel slides on the horizontal and vertical axes along the image, so for each pixel m × n multiplications are required; therefore the whole 2D filtering would have approximately a complexity of O(MNmn).It is easy to calculate that the total number of additions are (N + n 2 )(M + m 2 ).
In the simpler case used to exemplify our implementation, the filter kernel has size 12 × 12, while the image is 21 × 21; thus for usual convolution, there will be 63,504 multiplications with 27,225 additions.In our approach, only 100 inner products are used, with 3 × 3 multiplications and 3 × 3 additions for each, that is 100 × 33 = 900 multiplications and 900 additions, plus 12 × 90 additions in the pre-processing stage and 7 × 10 additions in the post-processing stage.As an additional example, for a larger value of the filter kernel where the filter kernel has the size 20 × 20 while the image is 35 × 35, 100 inner products are used, with 5 × 5 multiplications and 5 × 5 additions for each inner product, that is 100 × 5 × 5 = 2500 multiplications and 2500 additions, plus 20 × 90 additions in pre-processing stage and 70 additions in the post-processing stage.

Conclusions
The proposed design method for 2D circular filter banks is entirely analytical, without using any global numerical optimization.The advantages of the proposed method compared to other works are: it yields a factored 2D frequency response; the designed CFB is parametric, with adjustable characteristics; the CFB components are solved easily for any choice of number of filters and selectivity; the proposed implementation using polyphase and block filtering leads to a low complexity filter structure.This approach solves the problem of designing an adjustable and efficient circular FB with imposed specifications.The obtained results prove that the proposed rotationally invariant filter banks can be used in decomposing a given image into its subband components.
Taking into account the simulation results on test images and the fact that the original image can be reconstructed very accurately at least from a visual, subjective point of view from its component images, the authors intend in future work to study and investigate whether such CFBs (either uniform or non-uniform) could be used in sub-band coding schemes.While practically and intuitively this would seem possible, the required mathematical conditions for perfect reconstruction will have to be investigated rigorously.Regarding the implementation part, the authors will also study how to choose the decimation factors for the input image and the filter kernel, in order to obtain a very efficient, optimal design, and to minimize the number of arithmetic operations.

Sensors 2023 ,Figure 3 .
Figure 3. 1D prototypes, characteristics and contour plots for the first eight components of the circular filter bank; (a) low-pass filter; (h) band-pass filters BP1BP7.

Figure 3 .Figure 4 .
Figure 3. 1D prototypes, characteristics and contour plots for the first eight components of the circular filter bank; (a) low-pass filter; (b-h) band-pass filters BP1-BP7.

Figure 5 .
Figure 5. Characteristics and contour plots of the components of the dyadic CFB; (a) LP filter; (b-d) BP filters; (e) LP filter.

Figure 4 .
Figure 4. 1D prototypes, characteristics and contour plots for the last three components of the circular filter bank; (a,b) band-pass filters BP8, BP9; (c) high-pass filter.

Figure 4 .
Figure 4. 1D prototypes, characteristics and contour plots for the last three components of the circular filter bank; (a,b) band-pass filters BP8, BP9; (c) high-pass filter.

Figure 5 .
Figure 5. Characteristics and contour plots of the components of the dyadic CFB; (a) LP filter; (b-d) BP filters; (e) LP filter.

Figure 5 .
Figure 5. Characteristics and contour plots of the components of the dyadic CFB; (a) LP filter; (b-d) BP filters; (e) LP filter.

Figure 7 .
Figure 7. Relative energies calculated for the 11 images resulting at the output of the uniform cular filter bank.

Figure 7 .
Figure 7. Relative energies calculated for the 11 images resulting at the output of the uniform circular filter bank.

Figure 7 .Figure 8 .
Figure 7. Relative energies calculated for the 11 images resulting at the output of the uniform circular filter bank.

Figure 8 .
Figure 8. Image analysis using the dyadic circular FB: (a) LP filtered; (b-d) BP filtered with BPF1, BPF2, BPF3, respectively; (e) HP filtered; (f,g) recovered image by summing the first two and three components, respectively; (h) recovered image by summing all five component.

Table 1 .
Relative sub-band energies (in %) for the 11 images resulting at the output of the uniform CFB.

Table 1 .
Relative sub-band energies (in %) for the 11 images resulting at the output of the uniform CFB.
Sensors 2023, 23, x FOR PEER REVIEW 13 of

Table 2 .
Relative sub-band energies (in %) for the five images resulting at the output of the dyadic CFB.

Table 2 .
Relative sub-band energies (in %) for the five images resulting at the output of the dyadic CFB.